#ifndef _CU_IMPROVED_PERLIN_NOISE_H_
#define _CU_IMPROVED_PERLIN_NOISE_H_

#include <cuda_runtime.h>
#include "../cuCommon/cuTypes.h"

__constant__ __device__ int   cu_p[512];
__constant__ __device__ float cu_g[48];


namespace cuImprovedPerlinNoise
{

    static int host_p[514] =  
    { 
        151,160,137,91,90,15,
        131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
        190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
        88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
        77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
        102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
        135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
        5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
        223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
        129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
        251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
        49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
        138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
        151,160,137,91,90,15,
        131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
        190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
        88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
        77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
        102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
        135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
        5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
        223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
        129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
        251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
        49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
        138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
    };

    static float host_g[] = 
    {
         1, 1, 0,
        -1, 1, 0,
         1,-1, 0,
        -1,-1, 0,  
         1, 0, 1,
        -1, 0, 1,
         1, 0,-1,
        -1, 0,-1,  
         0, 1, 1,
         0,-1, 1,
         0, 1,-1,
         0,-1,-1,
         1, 1, 0,
         0,-1, 1,
        -1, 1, 0,
         0,-1,-1,
    };

    SD_CALL float  fade(float t);
    SD_CALL float  lerp(float t, float a, float b);
    SD_CALL float  dot(float a, float b, float c, 
                       float x, float y, float z);

    SD_CALL float  cuBrownianMotion   ( float3 point, 
                                        int    octaves    = 4, 
                                        float  lacunarity = 2.0, 
                                        float  gain       = 0.5 );

    SD_CALL float3 cuBrownianMotionf3 ( float3 point, 
                                        int    octaves    = 4, 
                                        float  lacunarity = 2.0f, 
                                        float  gain       = 0.5f);

    SD_CALL float  cuBrownianMotion   ( float4 point, 
                                        int    octaves    = 4, 
                                        float  lacunarity = 2.0, 
                                        float  gain       = 0.5 );

    SD_CALL float3 cuBrownianMotionf3 ( float4 point, 
                                        int    octaves    = 4, 
                                        float  lacunarity = 2.0f, 
                                        float  gain       = 0.5f);
    SD_CALL float  cuTurbulence       ( float3 point, 
                                        int    octaves    = 4, 
                                        float  lacunarity = 2.0f, 
                                        float  gain       = 0.5f);

    SD_CALL float cuTurbulence        ( float4 point, 
                                        int    octaves    = 4, 
                                        float  lacunarity = 2.0f, 
                                        float  gain       = 0.5f);
    SH_CALL  
    void   cuInit()
    {
        cudaMemcpyToSymbol("cu_p", host_p, 512 * sizeof( int   ));
        cudaMemcpyToSymbol("cu_g", host_g,  48 * sizeof( float ));
    }

    SD_CALL 
    float fade(float t) 
    { 
        return t * t * t * (t * (t * 6.0f - 15.0f) + 10.0f); 
    }

    SD_CALL 
    float lerp(float t, float a, float b) 
    { 
        return a + t * (b - a); 
    }

    SD_CALL 
    float dot(float a, float b, float c, 
              float x, float y, float z) 
    { 
        return a*x + b*y + c*z; 
    }

    SD_CALL 
    float cuGrad(int   hash, 
                 float x, 
                 float y,
                 float z) 
    {
        int h = (hash & 15) * 3;
        return dot(cu_g[h+0], cu_g[h+1], cu_g[h+2], x, y, z);
    }

    SD_CALL 
    float cuGrad(int   hash, 
                 float x, 
                 float y, 
                 float z, 
                 float w) 
    {
        int h = hash & 31;                                           // CONVERT LO 5 BITS OF HASH 
                                                                     // TO 32 GRAD DIRECTIONS.
        float a=y, b=z, c=w;                                         // X,Y,Z

        switch (h >> 3) 
        {                                                            // OR, DEPENDING ON HIGH ORDER 2 BITS:
            case 1: a=w;b=x;c=y;break;                               // W,X,Y
            case 2: a=z;b=w;c=x;break;                               // Z,W,X
            case 3: a=y;b=z;c=w;break;                               // Y,Z,W
        }

        return ((h&4)==0 ? -a:a) + ((h&2)==0 ? -b:b) + ((h&1)==0 ? -c:c);
    }

    SD_CALL 
    float cuNoise(float x, 
                  float y, 
                  float z, 
                  float w)
    {
        int X = (int)floor(x) & 255;                                 // FIND UNIT HYPERCUBE
        int Y = (int)floor(y) & 255;                                 // THAT CONTAINS POINT.
        int Z = (int)floor(z) & 255;
        int W = (int)floor(w) & 255;

        x -= floor(x);                                               // FIND RELATIVE X,Y,Z,W
        y -= floor(y);                                               // OF POINT IN CUBE.
        z -= floor(z);
        w -= floor(w);

        float a = fade(x);                                           // COMPUTE FADE CURVES
        float b = fade(y);                                           // FOR EACH OF X,Y,Z,W.
        float c = fade(z);
        float d = fade(w);
        int A = cu_p[X  ]+Y, AA = cu_p[A]+Z, AB = cu_p[A+1]+Z;       // HASH COORDINATES OF
        int B = cu_p[X+1]+Y, BA = cu_p[B]+Z, BB = cu_p[B+1]+Z;       // THE 16 CORNERS OF
        int AAA = cu_p[AA]+W, AAB = cu_p[AA+1]+W;                    // THE HYPERCUBE.
        int ABA = cu_p[AB]+W, ABB = cu_p[AB+1]+W;
        int BAA = cu_p[BA]+W, BAB = cu_p[BA+1]+W;
        int BBA = cu_p[BB]+W, BBB = cu_p[BB+1]+W;
                                                                     // INTERPOLATE DOWN.
        return lerp(d, lerp(c, lerp(b,lerp(a,cuGrad(cu_p[AAA  ], x  , y  , z  , w  ), 
                                             cuGrad(cu_p[BAA  ], x-1, y  , z  , w  )),
                                      lerp(a,cuGrad(cu_p[ABA  ], x  , y-1, z  , w  ), 
                                             cuGrad(cu_p[BBA  ], x-1, y-1, z  , w  ))),
                               lerp(b,lerp(a,cuGrad(cu_p[AAB  ], x  , y  , z-1, w  ), 
                                             cuGrad(cu_p[BAB  ], x-1, y  , z-1, w  )),
                                      lerp(a,cuGrad(cu_p[ABB  ], x  , y-1, z-1, w  ),
                                             cuGrad(cu_p[BBB  ], x-1, y-1, z-1, w  )))),
                        lerp(c,lerp(b,lerp(a,cuGrad(cu_p[AAA+1], x  , y  , z  , w-1), 
                                             cuGrad(cu_p[BAA+1], x-1, y  , z  , w-1)),
                                      lerp(a,cuGrad(cu_p[ABA+1], x  , y-1, z  , w-1), 
                                             cuGrad(cu_p[BBA+1], x-1, y-1, z  , w-1))),
                               lerp(b,lerp(a,cuGrad(cu_p[AAB+1], x  , y  , z-1, w-1), 
                                             cuGrad(cu_p[BAB+1], x-1, y  , z-1, w-1)),
                                      lerp(a,cuGrad(cu_p[ABB+1], x  , y-1, z-1, w-1),
                                             cuGrad(cu_p[BBB+1], x-1, y-1, z-1, w-1)))));
    }

    SD_CALL 
    float cuNoise(float x, 
                  float y,
                  float z) 
    {
        int X = (int)floor(x) & 255;                                 // FIND UNIT CUBE THAT
        int Y = (int)floor(y) & 255;                                 // CONTAINS POINT.
        int Z = (int)floor(z) & 255;
        x -= floor(x);                                               // FIND RELATIVE X,Y,Z
        y -= floor(y);                                               // OF POINT IN CUBE.
        z -= floor(z);

        float u = fade(x);                                           // COMPUTE FADE CURVES
        float v = fade(y);                                           // FOR EACH OF X,Y,Z.
        float w = fade(z);

        int A = cu_p[X  ]+Y, AA = cu_p[A]+Z, AB = cu_p[A+1]+Z;       // HASH COORDINATES OF
        int B = cu_p[X+1]+Y, BA = cu_p[B]+Z, BB = cu_p[B+1]+Z;       // THE 8 CUBE CORNERS,

        return lerp(w, lerp(v, lerp(u, cuGrad(cu_p[AA  ], x  , y  , z   ),   // AND ADD
                                       cuGrad(cu_p[BA  ], x-1, y  , z   )),  // BLENDED
                               lerp(u, cuGrad(cu_p[AB  ], x  , y-1, z   ),   // RESULTS
                                       cuGrad(cu_p[BB  ], x-1, y-1, z   ))), // FROM  8
                       lerp(v, lerp(u, cuGrad(cu_p[AA+1], x  , y  , z-1 ),   // CORNERS
                                       cuGrad(cu_p[BA+1], x-1, y  , z-1 )),  // OF CUBE
                               lerp(u, cuGrad(cu_p[AB+1], x  , y-1, z-1 ),
                                       cuGrad(cu_p[BB+1], x-1, y-1, z-1 ))));
    }

    SD_CALL float cuNoise(float3 p)
    {
        return cuNoise(p.x, p.y, p.z);
    }

    SD_CALL float cuNoise(float4 p)
    {
        return cuNoise(p.x, p.y, p.z, p.w);
    }

    // vector Noise
    SD_CALL float3 cuNoise3f(float3 p)
    {
        return f3(cuNoise(p), 
                         cuNoise(p + f3(127, 71,  59)), 
                         cuNoise(p + f3( 37, 11, 131)));
    }

    SD_CALL float3 cuNoise3f(float4 p)
    {
        return f3(cuNoise(p), 
                         cuNoise(p + f4(127, 71,  59, 19)), 
                         cuNoise(p + f4( 37, 11, 131, 61)));
    }

    // Brownian Motion
    SD_CALL float cuBrownianMotion( float3 p, 
                                    int    octaves, 
                                    float  lacunarity, 
                                    float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;    
        for(int i=0; i<octaves; i++) 
        {
            noiseSum += cuNoise(p * freq) * amp;
            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }

    SD_CALL 
    float3 cuBrownianMotionf3( float3 p, 
                               int    octaves, 
                               float  lacunarity, 
                               float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float3 noiseSum = {0.0f, 0.0f, 0.0f};    
        for(int i=0; i<octaves; i++) 
        {
            noiseSum += cuNoise3f(p * freq) * amp;
            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }

    SD_CALL 
    float cuTurbulence( float3 p, 
                        int    octaves, 
                        float  lacunarity, 
                        float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;
        for(int i=0; i<octaves; i++) 
        {
            noiseSum += abs(cuNoise(p * freq))*amp;

            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }

    // Brownian Motion
    SD_CALL 
    float cuBrownianMotion( float4 p, 
                            int    octaves, 
                            float  lacunarity, 
                            float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum  = 0.0f;    
        for(int i=0; i<octaves; i++) 
        {
            noiseSum  += cuNoise(p * freq) * amp;

            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }

    SD_CALL float3 cuBrownianMotionf3( float4 p, 
                                       int    octaves, 
                                       float  lacunarity, 
                                       float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float3 noiseSum = {0.0f, 0.0f, 0.0f};    
        for(int i=0; i<octaves; i++) 
        {
            noiseSum += cuNoise3f(p * freq) * amp;
            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }

    SD_CALL float cuTurbulence( float4 p, 
                                int    octaves, 
                                float  lacunarity, 
                                float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;
        for(int i=0; i<octaves; i++) 
        {
            noiseSum += abs(cuNoise(p * freq)) * amp;

            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }
};

#endif